/** ANAYLSIS_12.C
 *  Analyse Macros zum F-Praktikumsversuch 09: Neutronennachweiss
 * 
 *  date:   22-Dec-2011
 *  author: Ch.Weinsheimer 
 */

#define N 1024

#include <iostream>



/**
 * User Defined Functions
 */

/// Define gauss function plus offset
double gaussian(double* vars, double* pars) {
    return (pars[3]+pars[0]*TMath::Gaus(vars[0],pars[1],pars[2]));
}

/// Define gauss function with linear background
double gaussian_back(double* vars, double* pars) {
    return (pars[0]*TMath::Gaus(vars[0],pars[1],pars[2]) + pars[3]*vars[0] + pars[4]);
}

/// Define three gauss function over linear background
double carbonfit(double* vars, double* pars) {
    return (pars[0]*TMath::Gaus(vars[0],pars[1],pars[2]) + pars[3]*TMath::Gaus(vars[0],pars[4],pars[5]) + pars[6]*TMath::Gaus(vars[0],pars[7],pars[8]) + pars[9]*vars[0] + pars[10]);
}

/// Define linear function
double lin(double* vars, double* pars) {
    return (pars[0]*vars[0]+pars[1]);
}




/// Read data from file and plot a histogram
void histo(char* filename)
{ 

    // File Opening, Reading and Closing
    FILE* file = fopen(filename, "r");
    Int_t channel[N], value[N];
    for(Int_t i=0; i < N; i++){
        fscanf(file, "%d %d", &channel[i],&value[i]);
    }
    fclose(file);


    // Create and fill 1-D Histogram
    TH1F* h1 = new TH1F("h1", "Data", N, 1, N);
    for(Int_t i=0; i < N; i++) h1->Fill(channel[i],value[i]);


    // Set styles
    gROOT->SetStyle("Plain");
    gStyle->SetOptStat(0);
    h1->SetTitle("Spectrum");
    h1->GetXaxis()->SetTitle("channel");
    h1->GetYaxis()->SetTitle("counts");

    // Create Canvas
    TCanvas* c1 = new TCanvas("c1", "Versuch09", 1024, 768);
    c1->cd(1);

    // Paint Graph
    h1->Draw("PL");
}



/// Read data from file and plot a histogram
void fitPeaks(char* filename,int const npeaks)
{

    // File Opening, Reading and Closing
    FILE* file = fopen(filename, "r");
    Int_t channel[N], value[N];
    for(Int_t i=0; i < N; i++){
        fscanf(file, "%d %d", &channel[i],&value[i]);
    }
    fclose(file);


    // Create and fill 1-D Histogram
    TH1F* h2 = new TH1F("h2", "Data_fitted", N, 1, N);
    for(Int_t i=0; i < N; i++) h2->Fill(channel[i],value[i]);


    // Set styles
    gROOT->SetStyle("Plain");
    gStyle->SetOptFit(111);
    gStyle->SetOptStat(0);
    h2->SetTitle("Fitted Spectrum");
    h2->GetXaxis()->SetTitle("channel");
    h2->GetYaxis()->SetTitle("counts");

    // Create Canvas
    TCanvas* c2 = new TCanvas("c2", "Versuch09", 1024, 768);
    c2->cd(1);



    // Declare Fit-Parameters and Fit-Funtions
    double mean[npeaks], stdv[npeaks], norm[npeaks], offset[npeaks], min[npeaks], max[npeaks];
    TF1 fit[npeaks];
    TH1F * hf[npeaks];  // Store Fitted Histograms here
    for (unsigned int j = 0; j < npeaks; j++) {
      cout << "Setting Parameters for Peak " << j + 1 << ": " << endl << "******************************" << endl;
      cout << "Position   : "; cin >> mean[j];
      cout << "Lower Bound: "; cin >> min[j];
      cout << "Upper Bound: "; cin >> max[j];
      
      char name[5];
      char namef[4];
      sprintf( name,"fit%d", j);
      sprintf( namef,"hf%d", j);

      hf[j] = (TH1F*)h2->Clone(namef);

      fit[j] = new TF1(name,gaussian_back,min[j],max[j],5);
      fit[j].SetParNames("Norm", "Mean", "StdDev", "m", "Offs");
      fit[j].SetParameters(100,mean[j],10,1,1);
      fit[j].SetLineColor(kRed);
      fit[j].SetLineWidth(2);

      hf[j]->Fit(&fit[j],"R","Sames");

      //h2->Fit(&fit[j],"R");
      cout << "Chi-Square: " << fit[j].GetChisquare() << "/" << fit[j].GetNDF() << endl;
    }

    // Update Canvas
    c2->Update();


    // Paint Graph
    h2->DrawClone("Same");
    for (j=0; j<npeaks; j++) {
      fit[j].DrawClone("Same");
    }
}



void energyCal()
{
    double meanCh[5] = {66.6, 82.4, 133.6, 146.3, 177.9};
    double energy[5] = {511.0, 622.0, 1173.0, 1274.0, 1333.0};
    double stdCh[5] = {0.9, 1.0, 1.6, 1.3, 36.0};

    // Create TGraph representation
    TGraphErrors* g1 = new TGraphErrors(5,meanCh,energy,stdCh,0);

    // Set styles
    gROOT->SetStyle("Plain");
    gStyle->SetOptFit(111);
    gStyle->SetOptStat(0);
    g1->SetTitle("Energy Calibration");
    g1->GetXaxis()->SetTitle("Channel");
    g1->GetYaxis()->SetTitle("Energy [keV]");

    // Create Canvas
    TCanvas* c3 = new TCanvas("c3", "Versuch09", 1024, 768);
    c3->cd(1);

    TF1 energyFit("energyFit", lin, 0, 150, 2);
    energyFit.SetParNames("Slope", "Offset");
    energyFit.SetParameters(10.0,1.0);
    energyFit.SetLineColor(kRed);
   
    g1->Fit(&energyFit,"0");

    
    energyFit->SetMarkerStyle(2);
    energyFit->SetMarkerSize(2);
    energyFit->SetLineWidth(2);
   
    g1->Draw("APE");
    energyFit.DrawClone("Same");
    
}


/// Plot Time Calibration Measurement
void timecal()
{
  Float_t peaks[16] = {36.0, 87.0, 136.0, 187.0, 231.0, 282.0, 331.0, 382.0, 425.0, 476.5, 526.0, 577.0, 621.0, 672.0, 721.0, 773.0};
  Float_t dpeaks[16] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
  Float_t times[16] = {4.0, 8.0, 12.0, 16.0, 20.0, 24.0, 28.0, 32.0, 36.0, 40.0, 44.0, 48.0, 52.0, 56.0, 60.0, 64.0};
  Float_t dtimes[16] = {0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5};

  // Create timegraph
  TGraphErrors * gr = new TGraphErrors(16,peaks,times,dpeaks,dtimes);
  gr->SetMarkerColor(kRed);
  gr->SetMarkerStyle(20);
  gr->SetMarkerSize(0.7);

  TF1 tFunc("lin", lin, 200, 500, 2);
  tFunc.SetParNames("Slope", "Offset");
  tFunc.SetParameters(1.0,0.0);
  tFunc.SetLineColor(kBlue);
  tFunc.SetLineWidth(2);
 
  gr->Fit(&tFunc);



  // Setze Styles
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);

  // Setze Beschriftungen und Styles
  gr->SetTitle("Time Calibration");
  gr->GetXaxis()->SetTitle("channel");
  gr->GetYaxis()->SetTitle("time / ns");

  // Erzeuge das Fenster
  TCanvas* Canva = gROOT->Find Object("Canva");
  if(!Canva) Canva = new TCanvas("Canva", "Versuch_12 TimeCal", 1024, 768);
  else Canva->Clear();
  Canva->cd(1); // Waehle das Feld an

  // Male Graph
  gr->Draw("APE");
  tFunc.Draw("Same");
}


/// Read data from file and plot a histogram
void carbon(char* filename)
{
    int const ub = 510;
    int const lb = 275;

    // File Opening Reading and Closing
    FILE* file = fopen(filename, "r");
    Int_t channel[N], value[N];
    for(Int_t i=0; i < N; i++){
        fscanf(file, "%d %d", &channel[i],&value[i]);
    }
    fclose(file);

    // Create 1-D Histogram
    TH1F* h3 = new TH1F("h3", "Fitted Histogramm", N, 250, 600);

    // Fill histogramm
    for(Int_t i=250; i < 600; i++){
      h3->Fill(channel[i],value[i]);
    }


    // Erzeuge das Fenster
    gROOT->SetStyle("Plain");
    TCanvas* c4 = new TCanvas("c4", "Fitted Versuch_12", 1024, 768);
    c4->cd(1); // Waehle das Feld an

    // Setze Styles
    gStyle->SetOptStat(0);
    gStyle->SetOptFit(111);


    // Declare Fit-Parameters and Fit-Funtions
    TF1 background("background","[0]*x+[1]",lb,ub);
    background.SetLineWidth(2);
    background.SetLineStyle(2);
    background.SetLineColor(kBlue);

    TF1 signal("signal","[0]*TMath::Gaus(x,[1],[2])+[3]*TMath::Gaus(x,[4],[5])+[6]*TMath::Gaus(x,[7],[8])",lb,ub);
    signal.SetLineWidth(2);
    signal.SetLineStyle(2);
    signal.SetLineColor(kRed);

    TF1 fit("fit",carbonfit,lb,ub,11);
    fit.SetLineWidth(2);
    fit.SetLineColor(kRed);
    double norm1 = 200.0; double mean1 = 363.4; double stddev1 = 15.0;
    double norm2 = 200.0; double mean2 = 415.3; double stddev2 = 15.0;
    double norm3 = 200.0; double mean3 = 475.0; double stddev3 = 15.0;
    double m = -2.0; double offset = 1000.0;
    fit.SetParameters(norm1,mean1,stddev1,norm2,mean2,stddev2,norm3,mean3,stddev3,m,offset);
    fit.SetParNames("Norm1","Mean1","StdDev1","Norm2","Mean2","StdDev2","Norm3","Mean3","StdDev3","m","Offset");

    // Perform Fit
    h3->Fit(&fit,"R");
    cout << "Chi-Square: " << fit.GetChisquare() << "/" << fit.GetNDF() << endl;


    background.SetParameters( fit.GetParameter(9), fit.GetParameter(10) );
    signal.SetParameters( fit.GetParameter(0), fit.GetParameter(1), fit.GetParameter(2),
			  fit.GetParameter(3), fit.GetParameter(4), fit.GetParameter(5),
			  fit.GetParameter(6), fit.GetParameter(7), fit.GetParameter(8)
			);

    // Setze Beschriftungen und Styles
    h3->SetTitle("Spectrum");
    h3->GetXaxis()->SetTitle("channel");
    h3->GetYaxis()->SetTitle("counts");
    h3->SetMarkerStyle(2);
    h3->SetLineStyle(0);
    
    //    hp->GetYaxis()->SetRangeUser(,1050);

    // Male Graph
    h3->Draw("PL");
    background.DrawClone("Same");
    signal.DrawClone("Same");


}



void timeOfFlight(char* filename, int const npeaks)
{
    // File opening, reading and closing
    FILE* file = fopen(filename, "r");
    Int_t channel[N], value[N];
    for(Int_t i=0; i < N; i++){
        fscanf(file, "%d %d", &channel[i],&value[i]);
    }
    fclose(file);

    // Create 1-D histogram
    TH1F* h5 = new TH1F("h5", "Fitted Histogramm", N, 0, N);

    // Fill histogramm
    for(Int_t i=0; i < N; i++){
      h5->Fill(channel[i],value[i]);
    }

    // Create the canvas for drawing
    gROOT->SetStyle("Plain");
    TCanvas* c5 = new TCanvas("c5", "Fitted Versuch_12", 1024, 768);
    c5->cd(1); // Waehle das Feld an

    // Set Styles
    gStyle->SetOptStat(0);
    gStyle->SetOptFit(111);

    // Background estimation
    TSpectrum * s = new TSpectrum(2);
    TH1 * h_back =  s->Background(h5);
    h_back->SetLineColor(kBlue);
    h_back->SetLineStyle(2);

    // Create new signal histogramm with the background subtracted
    TH1F * signal = (TH1F*)h5->Clone("signal");
    signal->Add(h_back,-1.0);
    signal->SetLineColor(kRed);

    // Perform fits with gaussian functions
    double mean[npeaks], stdv[npeaks], norm[npeaks], offset[npeaks], min[npeaks], max[npeaks];
    TF1 fit[npeaks];
    TH1F * hf[npeaks];  // Store Fitted Histograms here
    for (unsigned int j = 0; j < npeaks; j++) {
      cout << "Setting Parameters for Peak " << j + 1 << ": " << endl << "******************************" << endl;
      cout << "Position   : "; cin >> mean[j];
      cout << "Lower Bound: "; cin >> min[j];
      cout << "Upper Bound: "; cin >> max[j];
      
      char name[5];
      char namef[4];
      sprintf( name,"fit%d", j);
      sprintf( namef,"hf%d", j);

      hf[j] = (TH1F*)signal->Clone(namef);

      fit[j] = new TF1(name,gaussian,min[j],max[j],4);
      fit[j].SetParNames("Norm", "Mean", "StdDev", "Offs");
      fit[j].SetParameters(100,mean[j],10,1);
      fit[j].SetLineColor(kBlue);
      fit[j].SetLineWidth(2);

      hf[j]->Fit(&fit[j],"R","Sames");

      //h2->Fit(&fit[j],"R");
      cout << "Chi-Square: " << fit[j].GetChisquare() << "/" << fit[j].GetNDF() << endl;
    }

    // Update Canvas
    c5->Update();


    // Draw
    h5->DrawClone("Same");
    for (j=0; j<npeaks; j++) {
      fit[j].DrawClone("Same");
    }
    signal->Draw("Same");
    h_back->Draw("Same");
}
